% Matlab code to visualize the equilibrium for the calibrated model

% To run this file seperately, specify the root directly below:
% cd 'C:/Users/Jilun/Dropbox/WorldEconomy_RiskIntolerance/paper RFS final submission/replication/numerical and calibration exercises'


% This code is for the case with LSAPs. It makes two figures
% FIGURE 7: Plot the equilibrium price for a range of productivites, without and with LSAPs
% FIGURE B2: Graphically illustrate the equilibrium, without and with LSAPs, for a particular productivity shock (included in the appendix) 


clear all;
% Baseline parameters (will change later to illustrate multiplicity)
global g rho sigma PStar taub tauh k l tau tauBarBench zStar
loadVariables_calibration; % Sets the above parameters and the baseline shocks

zHigh = 1;
% Find the threshold that eliminates the good equilibrium
zLow = zStar-0.03; 
zPost = 0.97;
lsap = 0.01;

ph = exp(rho+g - tauh/sigma^2);
 
% Recall that zh=l/ph
    %zh = l/ph;
    %zLow= zh;

zArr = [zLow:0.0001:zHigh];
pArr = zArr*0;
% Solve the equilibrium for a range of z's
for i=1:length(zArr)
    z0= zArr(i);
    eq{i} = solveEquilibrium(z0,0);
    eq_lsap{i} = solveEquilibrium(z0,lsap);
    pArr(i) = eq{i}.p;
    pArr_lsap(i) = eq_lsap{i}.p;
    uniqueArr(i) = eq{i}.unique; % Whether the equilibrium is unique
    uniqueArr_lsap(i) = eq_lsap{i}.unique; % Whether the equilibrium is unique
    % When there are multiple equilibria, the code automatically picks the
    % highest equilibrium
end

% With these parameters, there is always a bad equilibrium
% Find the first z that also allows for a good equilibrium
zInd= find(uniqueArr==0,1,'first'); 
zInd_lsap= find(uniqueArr_lsap==0,1,'first'); 


%%%%%%%%%%%% FIGURE 7
%%%%%%%%%%%% Calibrated equilibrium and its approximation without LSAPs

figure('Position',[30 30 830 430]);
hold on;
% Use the Matlab default blue color, which is [0 0.4470 0.7410]
h(1) = plot(zArr(zInd_lsap:end),pArr_lsap(zInd_lsap:end),'LineStyle','-','LineWidth',3);
defaultBlue = get(h(1),'Color');
h(2) = plot(zArr(zInd:end),pArr(zInd:end),'LineStyle','--','Color',defaultBlue,'LineWidth',3);

xlabel('Productivity, $z_0$','Interpreter','Latex','FontSize',20);
title('Price per productivity, $p_0$ (in the best equilibrium)','Interpreter','Latex','FontSize',20);
zLow = zArr(zInd_lsap); zHigh =zStar+0.005;
pLow = pArr(zInd); pHigh = pArr(end);
xlim([zLow zHigh]); ylim([pLow pHigh]);
set(gca,'Xtick',[0:0.005:1]);

% Plot a line for the equilibrium shock discussed in the main text
indPost = find(zArr>=zPost,1,'first');
pPost = pArr(indPost);
pPost_lsap = pArr_lsap(indPost);
pChange = (pPost_lsap - pPost)/lsap;
line([zPost zPost],[pLow pPost_lsap]);

plot(zPost,pPost,'ko-','LineWidth',3);
plot(zPost,pPost_lsap,'ko-','LineWidth',3);

legendStr1 = sprintf('With LSAP, $\\lambda=%1.2f$',lsap);
legendStr2 = sprintf('Without LSAP');
hL = legend(h([1 2]),legendStr1,legendStr2);
set(hL,'orientation','vertical','fontsize',15,'interpreter','latex');

saveas(gcf,[pwd '/figures/fig7.eps'],'epsc');


%%%%%%%%%%%% FIGURE B2
%%%%%%%%%%%% Actual and required Sharpe ratios for the calibrated equilibrium
%%%%%%%%%%%% illustrate a productivity shock that induces a demand recession, without and with LSAPs
eq = solveEquilibrium(zPost,0,1000); % No LSAP
eq_g = solveEquilibrium(zPost,lsap,1000);

pLow = 0.9;
% Plot for pArr>pLow

pArr = eq.pArr;
sharpeActualArr = eq.sharpeActualArr;
sharpeRequiredArr = eq.sharpeRequiredArr;
ind = eq.ind;

pArr_g = eq_g.pArr;
sharpeActualArr_g = eq_g.sharpeActualArr;
sharpeRequiredArr_g = eq_g.sharpeRequiredArr;
ind_g = eq_g.ind;
ph_g = eq_g.ph;
% Now plot the equilibrium curves

figure('Position',[30 30 830 430]);
hold on;
indRange = find(pArr>=pLow & pArr<=1);
h(1) = plot(pArr(indRange),sharpeActualArr(indRange),'b','LineWidth',3);
h(2) = plot(pArr(indRange),sharpeRequiredArr(indRange),'r--','LineWidth',3);
% Find and mark the equilibrium
plot(pArr(ind),(sharpeActualArr(ind)+sharpeRequiredArr(ind))/2,'ko-','LineWidth',3);

indRange_g = find(pArr_g>=pLow & pArr<=1);
h(3) = plot(pArr_g(indRange_g),sharpeRequiredArr_g(indRange_g),'r','LineWidth',3);
% Find and mark the equilibrium
plot(pArr_g(ind_g),(sharpeActualArr_g(ind_g)+sharpeRequiredArr_g(ind_g))/2,'ko-','LineWidth',3);

titleStr = sprintf('Large-scale asset purchases in a calibrated equilibrium');
title(titleStr,'Interpreter','Latex','FontSize',20);
xlabel('Normalized price, $p_0$','Interpreter','Latex','FontSize',20);

legendStr1 = sprintf('Actual Sharpe ratio, $r^f_0=0$');
legendStr2 = sprintf('Required Sharpe ratio without LSAPs');
legendStr3 = sprintf('Required Sharpe ratio with LSAPs, $\\lambda=%1.2f$',lsap);
hL = legend(h([1 3 2]),legendStr1,legendStr3,legendStr2);
set(hL,'orientation','vertical','fontsize',15,'interpreter','latex');
xlim([pLow 1]);
saveas(gcf,[pwd '/figures/figB2.eps'],'epsc');
